\documentclass{report}
\usepackage[UTF8]{ctex}
\usepackage{amsmath}
\usepackage{CJKfntef,pifont}
\usepackage[a4paper,left=3.14cm,right=2.5cm,top=2.5cm,bottom=2.5cm]{geometry}%页面设置宏包
\usepackage{graphicx,caption,url,float,bm,hyperref,mathrsfs,fancyhdr,listings,extarrows}
\hypersetup{colorlinks=false,
	 pdfborder={0 0 1}}
\usepackage[dvipsnames]{xcolor}
% 使用水印
%\usepackage{draftwatermark}
%\SetWatermarkText{内部资料}
%\SetWatermarkFontSize{3cm}% 最大5cm
% 设置页眉
\pagestyle{fancy}
%\lhead{\includegraphics[scale=0.08]{SchoolBadge2.jpg}}
%\lhead{湖南农大$ \bullet $暑期学校\qquad 内部资料$ \bullet $请勿外传}

% 控制目录的深度为1
%\setcounter{tocdepth}{1}
%\setcounter{secnumdepth}{1}

\captionsetup{figurename=图}

%代码风格设置
\lstset{language = R, basicstyle = {\ttfamily\small},keywordstyle = \color{RoyalBlue},frame = lines, commentstyle=\color{gray}}
%\lstset{numbers=none,
%	language=R,
%	numberstyle=\tiny,
%	keywordstyle=\color{blue!70}, commentstyle=\color{red!50!green!50!blue!50},
%	% frame=shadowbox,
%	rulesepcolor=\color{red!20!green!20!blue!20},
%	basicstyle=\ttfamily\small,
%	escapeinside=`', %将中文放在`和'之间，可以在lstlisting环境中准确的显示中文
%}

%opening
\title{R语言讲义}
\author{陈普}
\date{2018.7}

\begin{document}
	\maketitle
	\tableofcontents

\include{BaseofR}
\include{BaseofEco}

\chapter{微观计量经济学}
\section{因果推断的理论框架}
\subsection{潜在结果的表述框架}
因果推断，必有干预，在干预下讨论因果。被干预的个体称为\textbf{处理组}，未被干预的个体称为\textbf{控制组}。

在个体被干预之前，有几个干预状态就有几个\textbf{干预结果}，干预实现后，只有一个潜在结果可以观测。

一些符号表述：譬如考察大学教育对个体收入的影响，干预就是大学教育，
\begin{itemize}
	\item 任意个体$ i $，存在两种干预状态，\textbf{\color{RoyalBlue}$ D_i=1 $表示完成了大学教育，$ D_i=0 $表示完成高中教育未完成大学教育。}
	\item 任何个体无论其是否真的完成大学教育，事前总有两个状态：完成或者不完成大学教育。每个状态对应一个潜在结果，\textbf{\color{RoyalBlue}$ Y_{1i} $表示$ D_i=1 $的潜在结果，$ Y_{0i} $表示$ D_i=0 $的潜在结果}；
	\item \textbf{\color{RoyalBlue}$ Y_i $表示事后个体的观测结果}，很显然，
	\begin{align*}
	Y_i&=D_iY_{1i}+(1-D_i)Y_{0i}\\
	&=\begin{cases}
	Y_{1i}\qquad \text{如果}D_i=1\\
	Y_0i\qquad \text{如果}D_i=0\
	\end{cases} 
	\end{align*}
	\item 因果效应就可以定义为两种潜在结果的比较，
	\begin{equation}\label{ME_1}
	 {\color{RoyalBlue}\tau_i=Y_{1i}-Y_{0i}}
	\end{equation}
	
	\begin{itemize}
		\item 因果效应是干预后同一时间、同一物理个体潜在结果的比较。
		\item 仅有一个个体无法得到个体因果效应。
		\item 因果效应的关键是把未观测到的潜在结果（即反事实结果）估计出来。
	\end{itemize}
\end{itemize}
\subsection{稳定性假设和分配机制}
稳定性假设：
\begin{itemize}
	\item 不同个体个体潜在结果不会有交互影响（如同群效应）；
	\item 干预水平对所有个体一致；
\end{itemize}

\vspace{2em}
分配机制：为什么有的人在干预组，有的人在控制组，怎么分配的？根据分配机制是否已知，可将分配机制分成两类：\textbf{\color{RoyalBlue}随机化实验}（分配机制由实验者控制）和\textbf{\color{RoyalBlue}观测研究}（分配机制未知）。

引入关键概念：\textbf{\color{RoyalBlue}非混淆性（条件独立性）}假设：
\[ (Y_{0i},Y_{1i})\perp D_i|X_i \]

也即潜在结果独立于分配机制，个体在干预组还是控制组与潜在结果无关。其中$ X_i $是干预前变量，他们不受干预变量影响，但往往决定干预，通常包含两类：一是个体属性，如性别；二是干预实施前就已经决定的变量，如培训前的收入水平。

\subsection{总体平均因果效应}
公式\eqref{ME_1}表达的是个体因果效应，我们往往感兴趣总体平均因果效应（Average Treatment Effect, ATE），
\[ \tau_{ATE}=E[Y_{1i}-Y_{0i}] \]

或者干预组平均因果效应（Average Treatment Effect for the Treated, ATT），
\[ \tau_{ATT}=E[Y_{1i}-Y_{0i}|D_i=1] \]

或者控制组平均因果效应（Average Treatment Effect for the Treated, ATC），
\[ \tau_{ATC}=E[Y_{1i}-Y_{0i}|D_i=0] \]

\textbf{虚拟变量回归能否得到因果效应？}如
\[ Y_i=\alpha+\tau D_i+\varepsilon_i \]

我们知道此时，
\[ \hat \tau^{ols}=\overline Y_t-\overline Y_c \xLongrightarrow{p} E[Y_i|D_i=1]-E[Y_i|D_i=0]=\tau^{ols}\]

其中$ \overline{Y}_t,\overline{Y}_c $分别是干预组和控制组的结果变量的均值。那么，
\begin{align*}
\tau^{ols}& =E[Y_i|D_i=1]-E[Y_i|D_i=0]\\
& = E[Y_{1i}|D_i=1]-E[Y_{0i}|D_i=0]\\
& = \underbrace{E[Y_{1i}-Y_{0i}|D_i=1]}_{\text{ATT}}+\underbrace{E[Y_{0i}|D_i=1]-E[Y_{0i}|D_i=0]}_{\text{选择性偏误}}
\end{align*}

选择性偏误：干预组和控制组的个体，他们在没有干预下的潜在结果是不一样的。类似地，有，
\begin{align*}
\tau^{ols}& = \underbrace{E[Y_{1i}-Y_{0i}|D_i=0]}_{\text{ATC}}+\underbrace{E[Y_{1i}|D_i=1]-E[Y_{1i}|D_i=0]}_{\text{选择性偏误}}\hspace{11em} \text{\color{RoyalBlue}ATC}\\
\tau^{ols}& = \underbrace{E[Y_{1i}-Y_{0i}}_{\text{ATE}}+\underbrace{E[Y_{0i}|D_i=1]-E[Y_{0i}|D_i=0]}_{\text{选择性偏误}}+\\
&(1-P[D_i=1])(\underbrace{E[Y_{1i}-Y_{0i}|D_i=1]-E[Y_{1i}|D_i=1]-E[Y_{0i}|D_i=0]}_{\text{两组因果效应差异}})\hspace{5em}\text{\color{RoyalBlue}ATE}
\end{align*}

因此，OLS结果要包含因果含义，必须要求选择性偏误为0。
\section{随机实验}
自然科学通过将其他因素控制住，仅让关心的变量变化，观测结果的变化以得到因果解释。但社会科学难以做到控制住其他因素不变，因此随机化实验应运而生，它不要求其他影响因素不变，只要求干预分配是随机化的，同样也可以获得因果效应。

为什么？因为随机化的关键作用在于\textbf{\color{RoyalBlue}可以平衡两组个体其他因素的分布是相同的，}从而两组个体具有可比性。或者从计量经济学的观点看，对于模型，
\[ Y_i=\beta_0+\beta_1D_i+\varepsilon_i \]

由于$ D_i $是随机的，自然与其他未观测因素$ \varepsilon_i $无关，因此就不存在内生性问题，$ \beta_1 $的估计是一致的。


随机化实验的缺陷：现实中，接受干预的个体往往是个人选择的结果，不是随机的结果，他们的平均因果效应不是总体平均因果效应。

\subsection{随机化实验的分类}
\begin{itemize}
	\item 伯努利实验：掷硬币形式决定个体的干预状态。可能会出现所有个体均在干预组或控制组的情况。为克服该缺陷，提出完全随机化实验。	
	\item 完全随机化实验：事先确定一个固定的干预个体数目，剩下的作为控制组。但可能某些个体特征对潜在结果有重要影响，为克服该缺陷，提出分层随机化。	
	\item 分层随机化：根据特征分层以后，再在层内实施完全随机化。此时如何计算总体平均因果效应？
	
	用$ N_t(j) $表示$ j $层接受干预个体的数量，$ N(j) $表示$ j $层个体总数，则$ \frac{N_t(j)}{N(j)} $表示$ j $层个体接受干预的概率，亦称倾向指数，则总体平均因果效应可表达为，
	\begin{align*}
	\tau_{ATE}=&E[Y_{1i}-Y_{0i}]\\
	=& E[E[Y_{1i}-Y_{0i}]|G_i]\qquad \text{全期望公式}\\
= &E[Y_{1i}-Y_{0i}|G_i=1]Pr(G_i=1)+E[Y_{1i}-Y_{0i}|G_i=2]Pr(G_i=2)+\cdots+\\
	&E[Y_{1i}-Y_{0i}|G_i=J]Pr(G_i=J) \qquad \text{期望公式展开}\\
 = 	&\sum_{j=1}^N\tau_{ATE}(j)p(j)
	\end{align*}
\end{itemize}


\section{匹配}
匹配方法也需要条件独立假设（CIA）。

匹配的思路在于\textbf{\color{RoyalBlue}对于干预组个体，在控制组中寻找特征相似的控制组个体与其匹配，从而用该控制组个体的结果作为干预组个体的反事实结果。}
\subsection{协变量匹配}
回忆干预组平均处理效应：
\[ \tau_{ATT}=E[Y_{1i}|D_i=1]-E[Y_{0i}|D_i=1] \]

其中$ E[Y_{0i}|D_i=1] $不可观测。但在\textbf{完全随机化实验}中，因$ (Y_{0i},Y_{1i})\perp D_i $，故有$ E[Y_{0i}|D_i=1]=E[Y_{0i}|D_i=0] $，所以可以用控制组观测结果$ E[Y_{i}|D_i=0] $作为干预组反事实结果$ E[Y_{0i}|D_i=1] $的估计。在\textbf{分层随机化实验}中，因$ (Y_{0i},Y_{1i})\perp D_i |X_i$，故有$E[Y_{0i}|X_i,D_i=1]=E[Y_{0i}|X_i,D_i=0]  $，那么，
\begin{align*}
\tau_{ATT} & =E[Y_{1i}-Y_{0i}|D_i=1]\\
& = E[E[Y_{1i}-Y_{0i}|X_i,D_i=1]|D_i=1]\qquad \text{全期望公式}\\
& = E[E[Y_{1i}|X_i,D_i=1]-E[Y_{0i}|X_i,D_i=1]|D_i=1]\qquad \text{展开}\\
& = E[E[Y_{1i}|X_i,D_i=1]-E[Y_{0i}|X_i,D_i=0]|D_i=1]\qquad \text{条件独立假设，CIA}\\
& = E[\underbrace{E[Y_{i}|X_i,D_i=1]-E[Y_{i}|X_i,D_i=0}_{\tau_X}]|D_i=1]\qquad \text{转化成可观测结果}
\end{align*}

其中$ \tau_X $就是相同协变量$ X_i $层内两组结果平均值之差。如果$ X $是离散变量，那么，
\[ \tau_{ATT}=E[\tau_X|D_i=1]=\sum_x\tau_xPr(X_i=x|D_i=1) \]

上述表明：首先按照特征进行层内匹配，估计$ \tau_X $，然后按$ X $在处理组中的分布$ Pr(X=x|D_i=1) $进行加权即可得到干预组的平均因果效应。实际上，不同因果参数，仅仅在于权重不同。如总体平均因果效应，ATT为
\[ \tau_{ATT}=E[\tau_X]=\sum_x\tau_xPr(X_i=x) \]


为了保证匹配可以实施，需要特征为$ X_i $的个体在干预组和控制组都存在，即，\textbf{\color{RoyalBlue}共同区间假设}，
\[ 0 <Pr(D_i=1|X_i)<1 \]

\subsection{倾向得分匹配}
\subsubsection{倾向指数、平衡指数和倾向指数定理}
当$ X_i $是高维时，很可能在控制组找不到对应于干预组的个体，共同区间假设难以满足。为解决该难题，提出倾向得分指数匹配。定义\textbf{\color{RoyalBlue}倾向指数为}，
\[ Pr(X_i)=Pr(D_i=1|X_i) \]

即具有特征$ X_i $的个体接受干预的可能性。随机化实验，该概率由研究者控制。观测研究中，该指数未知，只能估计。

定义\textbf{\color{RoyalBlue}平衡指数}：

如果$ b(x) $满足$ D_i\perp X_i|b(X_i) $，则$ b(X_i) $是一个平衡指数。

通俗地说，基于平衡指数，分配与个人特征无关，也即干预组和控制组自变量分布相同。倾向指数就是一个平衡指数。实际上，你可以对你估计的倾向指数进行平衡指数的特征检验，估计的倾向指数没问题，应该满足协变量分布相同的特征。


\textbf{\color{RoyalBlue}倾向指数定理：如果有$ (Y_{0i},Y_{1i})\perp D_i|X_i $，则有，$ (Y_{0i},Y_{1i})\perp D_i|Pr(X_i) $}

\subsubsection{倾向指数匹配步骤：}
第一，定义相似性；第二，实施匹配；第三，评价匹配的效果；第四，估计因果效应。前两个步骤，需要反复进行，已达到最好的效果。
\paragraph{1. 定义相似性}
	\begin{itemize}
		\item 选哪些自变量？主要依据CIA假设，控制了变量，干预组与控制组个体无未观测差异。尽量多的引入，最终考察其引入后是否满足平衡指数特征。
		\item 如何测度相似性？
		\begin{itemize}
			\item 欧式距离？马氏距离？无论哪个距离，注意数据的标准化。
			\item 或直接利用倾向指数定义，即$ d(p(X_i),p(X_j)) $，但实践中往往采用线性化的倾向指数，即
			\[ l_i=\ln\frac{p(X_i)}{1-p(X_i)} \]
			否则，在倾向指数为0.3和0.7时，变化相同的单位，协变量的变化在0.7时会更大。线性化的倾向指数则无此问题，此时可以计算距离$ d(l_i,l_j) $。而倾向指数的估计，通常采用logit模型。
			\item 实践中选哪个，取决于手边的数据：协变量离散，可以精确匹配。协变量连续，只能使用（线性化）的倾向指数进行不精确匹配，它需要确定一个\textbf{卡尺}，只有在此卡尺范围内，个体才是相似的。
			
		\end{itemize}
	\end{itemize}
\paragraph{2. 实施匹配}
\begin{itemize}
	\item 近邻匹配：一对一（又称最近邻匹配），一对多（又称$ k $近邻匹配，数据多时考虑它）。注意：
	\begin{itemize}
		\item 如果存在距离相同多个控制个体怎么办？要么随机选。要么根据倾向指数从高到低进行排序，然后匹配，因为倾向指数高的难匹配（干预组倾向指数高）。要么根据多个控制个体的平均值作为干预组个体的匹配。
		\item 一对一是总要找一个个体与干预组匹配，但很可能此时差异太大，估计偏差加大。此时事先设定卡尺，卡尺范围内寻找匹配，最终可能会丢弃部分干预组样本。
		\item 控制组个体是否可以重复匹配？样本容量较小，则可考虑。因为重复匹配会降低估计精度。
		\item 最近邻匹配保证每一对距离最近，但对全部干预组个体而言，并非总体上最近。因此，有一种\textbf{\color{RoyalBlue}最优匹配}，总体上对所有干预组个体同时匹配，寻找总体距离最小。但计算耗时，对于平均因果效应估计收益不是太大。
		\item 半径匹配。类似于卡尺匹配，在半径内定义为相似，否则不相似。
	\end{itemize}
	\item 分层匹配：根据协变量或者倾向指数进行分层，Rosenbaum and Rubin (1985)表明根据倾向指数分成5个区间，减少偏差90\%。那么分几层是合适的呢？Imbens and Rubin (2015)提出了一种$ t $检验方法，检验层内的线性化指数是否无差异。
	
	具体地，令$ j $层干预组个体个数为$ N_{tj} $，控制组个体个数为$ N_{cj} $，则计算$ j $层内线性化倾向指数平均值为$ \bar l_{tj}=(\sum_{i:D_i=1}\hat l_i)/N_{tj},\bar l_{cj}=(\sum_{i:D_i=0}\hat l_i)/N_{cj} $，然后构造$ t $统计量，
	\[ t_j=\frac{\bar l_{tj}-\bar l_{cj}}{\sqrt{s_{l_j}^2(1/N_{tj}+1/N_{cj})}} \]
	
	若$ t $值大于1，则不够相似，根据层内中位数进一步分成新层，直到每层的$ t $都小于1。	
\end{itemize}

\paragraph{3. 评价匹配的效果}
\begin{itemize}
	\item 标准化平均值差异。
	\[ \hat\Delta_{ct}=\frac{\bar X_t-\bar X_c}{\sqrt{(s_t^2+s_c^2)/2}} \]
	
	其中，$ \bar X_t,\bar X_c $分别为某协变量的干预组和控制组平均值，$ s_t^2,s_c^2 $表示某协变量的样本方差。可见差异越小，该值越倾向于0。
	\item 对数标准差比。$ \Delta_{ct} $度量了一阶矩差异，对数标准差比度量的是二阶矩的差异，定义为，
	\[ \hat \varGamma_{ct}=\ln s_t-\ln s_c \]
	越接近0越好。
	\item 倾向指数的标准化平均值差异。非正态分布，前面二阶矩也不能代表整个分布，Imbens and Rubin (2015) 证明若两组倾向指数的期望值相同，则两组个体的协变量分布相同。那么就可以构造倾向指数的标注化平均值差异，
	\[ \hat\Delta_{ct}=\frac{\bar l_t-\bar l_c}{\sqrt{(s_{l,t}^2+s_{l,c}^2)/2}} \]
		
	\item 图形方法：
	\begin{itemize}
		\item 倾向指数分布图：画出干预组和控制组的倾向指数分布图。
		\item 分位数分布图：干预组和控制组按分位由低到高分别画在横轴和纵轴上，若相似，则是45度线。
		\item 标准化均值差异变化图：每个协变量匹配前后均值和标准差差异以图形呈现。
	\end{itemize}
\end{itemize}

\paragraph{4. 估计因果效应}
\textbf{匹配估计量}（除分层匹配估计量外）总可以写为，
\[ \hat \tau_{ATT}=\frac{1}{N_t}\sum_{i:D_i=1}\left[Y_i-\sum_{j\in M_j(i)}w(i,j)Y_j\right] \]

其中$ M_j(i) $表示对应于干预组个体$ i $的控制组个体集合。不同的匹配方法区别在与权重$ w(i,j) $不同。

\begin{itemize}
	\item 一对一匹配：$ w(i,j)=1 $。
	\item 一对$ k $匹配：$ w(i,j)=\frac{1}{k} $。
	\item 半径匹配：$ w(i,j)=\frac{1}{\#M_{j(i)}} $,其中$ \#M_{j(i)} $表示集合$ M_{j(i)} $内的个体数目。
	\item 核匹配：与干预组个体越近，权重越大。令$ K(X_i-X_j) $为核函数，则
	\[ w(i,j)=\frac{K(X_i-X_j)}{\sum_{j=1}^{\#M_{j(i)}}K(X_j-X_i)} \]
\end{itemize}

除非精确匹配，各种非精确匹配，其个体总是存在一定差异，这种差异是可以调整的，经过调整的比没经过调整的要好（Abadie and Imbens, 2011）。可利用匹配样本如下估计(没有匹配上的控制组和处理组个体均不纳入)，
\[ Y_i=\alpha_c+\tau D_i+X'_i\beta_c+v_i \]

其中$ \beta_c $就是匹配偏误，$ \tau $就是经调整的匹配估计量。

\textbf{分层匹配估计量}：
\begin{itemize}
	\item 第一步，譬如第$ j $层，在层内对控制组和干预组分别取平均，然后再对这些平均值差分，得到该层的平均因果效应估计，
	\[ \hat\tau^{dif}(j)=\overline{Y}_t(j)-\overline{Y}_c(j)\qquad j=1,2,\cdots,J \]
	\item 第二步，总体平均因果效应：各层样本占比作权重，作一个各层平均因果效应的加权和，
	\[ \hat\tau_{ATE}^{strat}=\sum_{j=1}^J\hat\tau^{dif}(j)\cdot \frac{N(j)}{N} \]
	干预组平均因果效应：干预组各层样本占比作权重，作一个各层平均因果效应的加权和，
		\[ \hat\tau_{ATT}^{strat}=\sum_{j=1}^J\hat\tau^{dif}(j)\cdot \frac{N_t(j)}{N_t} \]
	\item 第三步，若愿意，也可以做一个偏误纠正，用$ j $层内所有样本做一个回归，
	\[ Y_i=\alpha(j)+\tau(j)D_i+X'_i\beta(j)+\varepsilon_i \]
	在对$ \hat\tau $做相应加权即可。
\end{itemize}

\textbf{逆概加权估计量}

观察到，
\begin{align*}
E\left[\frac{D_iY_i}{p(X_i)}\right]&=E\left\{E\left[\frac{D_iY_i}{p(X_i)}|X_i\right]\right\}\qquad\text{全期望公式}\\
	& = E\left\{\frac{1}{p(X_i)}E\left[D_iY_i|X_i\right]\right\}\qquad X_i\text{为条件，则}p(X_i)\text{看作常数}\\
	& = E\left\{\frac{1}{p(X_i)}E\left[Y_i|X_i,D_i=1\right]p(X_i)\right\}\qquad \text{全期望公式}E\left[D_iY_i|X_i\right]=E[E\left[D_iY_i|X_i,D_i\right]]\\
	& = E\{E[Y_i|D_i=1,X_i]\} = E\{E[Y_{1i}|X_i]\} \qquad\text{换成潜在结果}\\
	& = E[Y_{1i}]\qquad \text{全期望公式}
\end{align*}
类似地，有，
\[ E\left[\frac{(1-D_i)Y_i}{1-p(X_i)}\right]=E[Y_{0i}] \]
则，总体平均因果效应为，
\begin{align*}
\tau_{ATE}=&E[Y_{1i}-Y_{0i}]\\
 =& E\left[\frac{D_iY_i}{p(X_i)}-\frac{(1-D_i)Y_i}{1-p(X_i)}\right]\\
 =& E\left[\frac{(D-p(X_i))Y_i}{p(X_i)(1-p(X_i))}\right]
\end{align*}
对应的样本等价物为，
\[ \hat\tau_{ATE}^{ipw}=\frac{1}{N}\sum_{i=1}^NY_i\frac{D_i-\hat p(X_i)}{\hat p(X_i)(1-\hat p(X_i))} \]

注意：逆概加权估计量若倾向指数存在误设，则会产生较大偏误，特别是在极端倾向指数附近（0或1），$ \tau_{ATE}^{ipw} $将趋向无穷大。

\subsubsection{匹配方法：条件独立性假设的检验}
条件独立性假设本身无法检验，但有一些必要条件还是可以检验的。
\begin{itemize}
	\item 伪结果方法：用一个事实上没有受到干预影响的变量作为伪结果，观察其因果效应是否为0。譬如使用干预实施前的一个结果变量。
	\item 伪干预方法：比较两个控制组的结果，观察其因果效应是否为0。
\end{itemize}

\subsubsection{实践}
\verb|PSW_1.R|

\section{DID}
\textbf{\color{RoyalBlue}适用场景：事前所有个体都没有受到政策干预，事后只有一组个体受到政策干预。}

\subsection{一个例子：垃圾焚化炉区位对住房价格的影响}

1978年开始传言要在北安德沃市兴建一座垃圾焚化炉,1981年开始动工。人们不禁想考察相对焚化炉的远近是否影响了周边房价。你可能会仅仅使用1981年的数据估计一个简单模型,
\[ rprice=\gamma_0+\gamma_1 nearine+u \]

其中$ nearine $表示在焚化炉3公里以内为1,否则为0。得到$ \hat \gamma_0=101307.5,\hat \gamma_1=-30688.27 $,此时截距就是远离焚化炉的住房均价，从系数估计看，越近房价就要越低。

上面的方程无法证明离焚化炉越近房价就越低,因为你如果使用1978年的数据做同样的回归,会得到,
\[ \widehat rprice=82517.23-18824.37 nearinc \]

这表明焚化炉本身就建立在了房价较低的地方。因此,你不能说因为建了焚化炉就降低了房价。关键在于1978至1981年$ nearinc$系数的变化。这两个系数之差为,
\[ \hat \delta_1=-30688.27-(-18824.37)=-11863.9 \]

这个估计值就是DID估计量。因为它可以表示成,
\[ \hat \delta_1=(\overline{rprice}_{81,nr}-\overline{rprice}_{81,fr} )-(\overline{rprice}_{78,nr}-\overline{rprice}_{78,fr} )\]

实际上,上述结果可以通过回归来得到，而且回归得到的结果还附带给出了标准误，
\[ rprice=\beta_0+\beta_1y81+\beta_2nearinc+\delta_1 y81\cdot nearinc+u \]

当然,我们也可以在回归中增加解释变量。

\textbf{回归得到DID的理解}

\begin{table}[H]
	\centering
	\begin{tabular}{cccc}
		\hline
		& Before($ y81=0 $) & After($ y81=1 $) & After-Before\\ \hline
		控制组($ nearinc=0 $) & $\beta_0$ & $\beta_0+\beta_1$ & $\beta_1$\\
		处理组($ nearinc=1 $) & $ \beta_0+\beta_2 $& $\beta_0+\beta_1+\beta_2+\delta_1$ &$ \beta_1+\delta_1$\\
		处理组-控制组 & $ \beta_2 $ & $ \beta_1+\delta_1 $ & $\delta_1$\\
		\hline
	\end{tabular}
\end{table}

\subsection{DID的假设条件}
也要满足如下的条件独立性：
\[ (Y_{1it},Y_{0it})\perp D_{it}|X_{it},U_{i},t \]

其中，$U_{i} $是不可观测时常变量。同时还有如下假设，
\begin{itemize}
	\item 关键假设，\textbf{\color{RoyalBlue}共同趋势假设}：干预组没接受干预，其结果的变动趋势与控制组的变动趋势相同，即
	\[ E[Y_{0i,t}-Y_{0i,t-1}|X_{it},D_i=1] = E[Y_{0i,t}-Y_{0i,t-1}|X_{it},D_i=0]\]
	或写为，
	\[ E[\Delta Y_{0i,t}|X_{it},D_i=1] = E[\Delta Y_{0i,t}|X_{it},D_i=0]\]
	评价：
	\begin{itemize}
		\item 回忆回归的结果要解释成因果效应，必须要求选择性偏误为0，即$ E[Y_{0i}|X_i,D_i=1]=E[Y_{0i}|X_i,D_i=0] $，差异在于水平变量还是差分变量。正是差分消除了$ U_i $。
		\item DID 类似于增量匹配，所以也可以用匹配估计法。
	\end{itemize}
\begin{figure}[H]
\includegraphics[scale=0.6]{DrawDID.png}
\caption{共同趋势假设}
\end{figure}
	\item 共同区间假设：
	\begin{align*}
	&Pr(D_i=1)>0\hspace{10em} \text{干预组和控制组都有观测值}\\
	&Pr(D_iT_i=1|X_{it}) < 1\hspace{7.5em} \text{有干预组个体也要有控制组个体，两期都要有值}
	\end{align*}
\item  外生性假设：即协变量不受政策干预影响。
\item 政策干预不能有交互影响或溢出效应。
\end{itemize}

\subsection{多期数据的DID：恐怖袭击对CBD经济的影响}
\subsubsection{问题}
9.11之后，人们是不是不敢在标志性建筑中办公？Abadie and Dermisi (2008)以芝加哥三座标志性建筑为中心，500以内的为阴影区域，是干预组。500米以外的是非阴影区域，为控制组。使用1996年2季度至2006年2季度的数据，进行了DID分析。两类建筑空置率走势图如下，图中可以看到干预组空置率在2003年会明显增大。
\begin{figure}[H]
\includegraphics[scale=0.6]{DID_AB1.png}
\caption{干预组和控制组平均空置率走势图}
\end{figure}
\subsubsection{估计}
主要的估计模型：
\[ VacancyRate_{it}=\tau\cdot (shadow_i\cdot post911_t)+f_t+\eta_i+\varepsilon_{it} \]

其中，$ VacancyRate_{it} $是$ i $楼$ t $期的空置率，$ shadow_i $在阴影区域内取1，$ post911_t $在2001年9月11日之后取，$ f_t,\eta_i $分别是时间和建筑的固定效应。为保证结果的稳健性，作者还估计了另外三个模型：
\begin{align*}
VacancyRate_{it}&=\tau\cdot (DistancetoAnchor_i\cdot post911_t)+f_t+\eta_i+\varepsilon_{it}\\
VacancyRate_{it}&=\tau\cdot (DistancetoNonShadowArea_i\cdot post911_t)+f_t+\eta_i+\varepsilon_{it}\\
VacancyRate_{it}&=\tau\cdot (height_i\cdot post911_t)+f_t+\eta_i+\varepsilon_{it}
\end{align*}

其中，$ DistancetoAnchor_i $表示到标志性建筑的距离，越远恐怖袭击可能性越小，空置率越小，预期为负。$ DistancetoNonShadowArea_i $表示阴影内建筑距离阴影边界的距离，越近恐怖袭击可能性越小，预期为正。$ height_i $表示建筑的高度，越高恐怖袭击概率越大，空置率越大，预期为正。注意，这是DID的拓展，它们不再是虚拟变量而是连续变量。
\begin{figure}[H]
	\includegraphics[scale=0.5]{DID_AB2.png}
	\caption{四个DID模型估计结果}
\end{figure}
\subsubsection{稳健性检验}

 \textbf{\color{RoyalBlue}1. 因为911的影响是持续的、累积的，不是偶发的}，因此重新估计如下模型：
\begin{align*}
VacancyRate_{it}&=\alpha\cdot (shadow_i\cdot post911_t)+\tau\cdot (shadow_i\cdot since911_t)+f_t+\eta_i+\varepsilon_{it}\\
VacancyRate_{it}&=\alpha\cdot (DistancetoAnchor_i\cdot post911_t)+\tau\cdot (DistancetoAnchor_i\cdot since911_t)+f_t+\eta_i+\varepsilon_{it}\\
VacancyRate_{it}&=\alpha\cdot (DistancetoNonShadowArea_i\cdot post911_t)+\tau\cdot (DistancetoNonShadowArea_i\cdot since911_t)+f_t+\eta_i+\varepsilon_{it}\\
VacancyRate_{it}&=\alpha\cdot (height_i\cdot post911_t)+\tau\cdot (height_i\cdot since911_t)+f_t+\eta_i+\varepsilon_{it}
\end{align*}

其中$ since911_t $是自911以来经过了星期的个数。实证结果$ \tau $的符号符合预期。

\textbf{\color{RoyalBlue}2. 共同趋势的检验}

只使用2001年以前的数据，然后以1998年为伪恐怖袭击时间点，重新估计上述方程，发现系数全部不显著。

\subsection{时变干预DID：多个个体干预时间不同}
定义一个新的干预变量（在多期中，仅有一个时期是1，回忆$ post911_t $，它是911之后全部是1），
\[ D_{it}=\begin{cases}
1\qquad \text{如果个体}i\text{在}t\text{期被处理}\\
0\qquad \text{如果个体}i\text{在}t\text{期没有被处理}\\
\end{cases} \]

那么，对于带有超前或者滞后一期$ D_{it} $的方程为，
\begin{equation}\label{DID1}
 Y_{it}=\gamma_i+\lambda_t+\beta_{-1}D_{it-1}+\beta_0D_{it}+\beta_{+1}D_{it+1}+\gamma X_{it}+\varepsilon_{it}
\end{equation}

根据\eqref{DID1}式，则$ t-1,t,t+1 $三个时期对$ Y $的预测为，
\begin{align*}
E[Y_{it-1}|w]&=\gamma_i+\lambda_{t-1}+\beta_{-1}D_{it-2}+\beta_0D_{it-1}+\beta_{+1}D_{it}+\gamma X_{it-1}+\varepsilon_{it-1}\\
E[Y_{it}|w]&=\gamma_i+\lambda_t+\beta_{-1}D_{it-1}+\beta_0D_{it}+\beta_{+1}D_{it+1}+\gamma X_{it}+\varepsilon_{it}\\
E[Y_{it+1}|w]&=\gamma_i+\lambda_{t+1}+\beta_{-1}D_{it}+\beta_0D_{it+1}+\beta_{+1}D_{it+2}+\gamma X_{it+1}+\varepsilon_{it+1}
\end{align*}

如何干预，其全部序列可以表达成，
\[ \{w^j\}=\{D_{it-1},D_{it},D_{it+1}\}=\begin{cases}
w^1=(0,0,0)\\
w^2=(1,0,0)\\
w^3=(0,1,0)\\
w^4=(0,0,1)
\end{cases} \]

在$w^3  $干预组序列下，则对三个时期$ Y $的预测为，
\begin{align*}
E[Y_{it-1}|w]&=\gamma_i+\lambda_{t-1}+\beta_{+1}+\gamma X_{it-1}\\
E[Y_{it}|w]&=\gamma_i+\lambda_t+\beta_0+\gamma X_{it}\\
E[Y_{it+1}|w]&=\gamma_i+\lambda_{t+1}+\beta_{-1}+\gamma X_{it+1}
\end{align*}
在$w^1  $控制组序列下，三个时期$ Y $的预测为，
\begin{align*}
E[Y_{it-1}|w]&=\gamma_i+\lambda_{t-1}+\gamma X_{it-1}\\
E[Y_{it}|w]&=\gamma_i+\lambda_t+\gamma X_{it}\\
E[Y_{it+1}|w]&=\gamma_i+\lambda_{t+1}+\gamma X_{it+1}
\end{align*}

可见：
\begin{itemize}
	\item $ \beta_{+1} $是$ t $期干预对$ t-1 $期的影响。$ \beta_{0} $是$ t $期干预对$ t $期的影响。$ \beta_{-1} $是$ t $期干预对$ t+1 $期的影响。
	\item 一般而言，$ \beta_{+1} $应该为0，除非存在预期效应。
\end{itemize}
\subsection{实践}
%两期数据的固定效应估计就是DID。
\verb|DID.R|

也可以使用stata程序，调用\verb|ssc install diff|安装\verb|diff|命令。或通过\verb|findit diff|搜寻然后安装。该命令内涵丰富，包括分位DID，PSM-DID。


\section{合成控制法}
	当前文献习惯使用比较研究来讨论政策干预。但是这至少存在两个问题：第一，比较组的选择存在某种任意性,譬如万一比较组和处理组并不十分相似呢？第二，若使用的是个体（非加总）数据，则传统的推断技术仅仅考虑了数据加总时存在的误差而没有考虑控制组再现反事实结果能力的误差。
\subsection{模型}
\subsubsection{符号表述}
\begin{table}[H]
	\begin{tabular}{lp{36em}}
		$ Y_{it}^N $：& 地区$ i\in \{1.\cdots,J+1\} $在时间$ t\in \{1,\cdots,T\} $无干扰下的结果变量，其中$T_0$是干扰前的一个时期，故$ 1\le T_0<T $；\\
		$ Y_{it}^I $:& 有干扰的结果变量，注意此时$ t\in\{T_0+1,\cdots,T\} $；\\
		$ \alpha_{it} $:& $ = Y_{it}^I-Y_{it}^N $；\\
		$ D_{it} $:&是否干预的示性变量；
	\end{tabular}
\end{table}
\subsubsection{模型设置}
\begin{equation}\label{1}
Y_{it}^N=\delta_t+\bm{\theta}_t\bm{Z}_i+\bm{\lambda}_t\bm{\mu}_i+\varepsilon_{it} 
\end{equation}


其中，$ \bm{Z}_i $是未被政策干扰时的可观测向量$(r\times 1)$，$ \bm{\theta}_t $是未知参数向量$ (1\times r) $，$ \bm{\lambda}_t $是不可观测共同因子$ (1\times F) $，$ \bm{\mu}_i $是因子载荷$ F\times 1 $。

若考虑$ J\times 1 $维的权重矩阵$ \bm{W}=(w_2,\cdots,w_{J+1}) $，其中，$ w_j\ge 0(j=2,\cdots,J+1),w_2+\cdots +w_{J+1}=1 $。那么合成控制的结果变量就为，
\[ \sum_{j=2}^{J+1}w_jY_{jt}=\delta_t+\bm{\theta}_t\sum_{j=2}^{J+1}w_j\bm{Z}_j +\bm{\lambda}_t\sum_{j=2}^{J+1}w_j\bm{\mu}_j+\sum_{j=2}^{J+1}w_j\varepsilon_{jt}
\]

假设若存在一个这样的向量$ (w_2^*,\cdots,w_{J+1}^*) $使得合成控制的结果可以完全模拟干扰组$ T_0 $之前的结果，即，
\begin{equation}\label{2}
\sum_{j=2}^{J+1}w_j^*Y_{j1}=Y_{11},\sum_{j=2}^{J+1}w_j^*Y_{j2}=Y_{12},\cdots,\sum_{j=2}^{J+1}w_j^*Y_{jT_0}=Y_{1T_0},\sum_{j=2}^{J+1}w_j^*\bm{Z}_{j}=\bm{Z}_{1}
\end{equation}


那么，附录B表明只要$ \sum_{t=1}^{T_0}\bm{\lambda}'_t\bm{\lambda}_t $非奇异，则，
\begin{equation}\label{synth3}
Y_{it}^N-\sum_{j=2}^{J+1}w_j^*Y_{jt}=\sum_{j=2}^{J+1}w_j^*\sum_{s=1}^{T_0}\bm{\lambda}_t\left(\sum_{n=1}^{T_0}\bm{\lambda}'_n\bm{\lambda}_n\right)^{-1}\bm{\lambda}'_s(\varepsilon_{js}-\varepsilon_{1s})-\sum_{j=2}^{J+1}w_j^*(\varepsilon_{jt}-\varepsilon_{1t})
\end{equation}

在标准条件下，只要干扰前的期数足够长，使得其相对于临时冲击相对较大，则\eqref{synth3}式的均值将接近0。\textbf{\color{RoyalBlue}由于\eqref{2}式在实践中不太可能精确保持},所以只能尽可能使得它保持。
\begin{itemize}
	\item 在每一个应用研究中，这种不匹配的大小是可以计算的，如果拟合太差，还是不要使用合成控制了。
	\item 另外，即便拟合得很好，如果简单线性模型在任何特定样本集中无法在整个区域(region)中成立（内插偏误较大），\textbf{\color{RoyalBlue}研究者最好把待选池限制在具有与干预地区具有类似特征的地区里面。}	
\end{itemize}

注意，\eqref{1}式实际上是DID模型的一般化，经典的DID模型$ \bm{\lambda}_t $是一个常数，所以时间上的差分，然后再个体差分就可以得到因果效应。

\subsubsection{具体算法}
令
\[  \bm{W}=(w_2,\cdots,w_{J+1})' \] 

其中，$ w_j\ge 0,(j=2,\cdots,J+1),w_2+\cdots+w_{J+1}=1 $。受干预地区的结果变量为$ Y_{1t} $，其他未干预地区为$ Y_{jt},j=2,\cdots,J+1 $。再令$ T_0\times 1$的向量$ \bm{K}=(k_1,\cdots,k_{T_0})' $定义一个线性组合，
\begin{equation}\label{4}
\bar Y_i^{\bm{K}}=\sum_{s=1}^{T_0}k_sY_{is} 
\end{equation}


譬如若$ k_1=\cdots=k_{T_0}=1/T_0 $，则$ \bar Y_i^{\bm{K}} $就代表了处理前时期的一个简单平均。继续考虑若存在$ M $个这样的$ \bm{K}_1\cdots,\bm{K}_M $定义了$ M $个这样的线性组合。此时，令
\[\bm{X}_1=(\bm{Z}'_1,\bar Y_1^{\bm{K}_1},\cdots,\bar Y_1^{\bm{K}_M})'  \]
是一个$ (k\times 1) $维（$ k=r+M $）的表示干预地区的处理前特征。类似地，未干预地区可以记为$ \bm{X}_0 $，其维度为$ (k\times J) $，具体地，
\[ \bm{X}_0=\begin{pmatrix}
\bm{Z}_2 & \bm{Z}_3 & \cdots & \bm{Z}_{J+1}\\
\bar Y_2^{\bm{K}_1} & \bar Y_3^{\bm{K}_1} & \cdots & \bar Y_{J+1}^{\bm{K}_1}\\
\vdots & \vdots &\ddots & \vdots\\
\bar Y_2^{\bm{K}_M} & \bar Y_3^{\bm{K}_M} & \cdots & \bar Y_{J+1}^{\bm{K}_M}
\end{pmatrix} \]

因此，可以选择一个向量$ \bm{W}^* $来最小化距离$ \parallel \bm{X}_1-\bm{X}_0\bm{W}\parallel $。那又如何确定$ \bar Y_i^{\bm{K}_1},\cdots,\bar Y_i^{\bm{K}_M} $？一个明显的选择就是令$ \bar Y_i^{\bm{K}_1}=Y_{i1},\cdots,\bar Y_i^{\bm{K}_{T_0}}=Y_{iT_0} $，实践中，为了计算上的简单，只需考虑处理前结果的某些线性组合是否可以令\eqref{2}式保持地最好。

另外，$ \parallel \bm{X}_1-\bm{X}_0\bm{W}\parallel $是欧式距离，也可以使用马氏距离，
\begin{equation}\label{5}
\parallel \bm{X}_1-\bm{X}_0\bm{W}\parallel_{\bm{V}}=\sqrt{(\bm{X}_1-\bm{X}_0\bm{W})'\bm{V}(\bm{X}_1-\bm{X}_0\bm{W})}
\end{equation}

其中$ \bm{V} $是一个$ (k\times k) $的对称正定矩阵。但要注意若因变量与$ \bm{X}_0,\bm{X}_1 $中的解释变量呈现高度非线性的关系，那么\textbf{\color{RoyalBlue}内插偏误可能会很大。此时目标函数应该是最小化$\parallel \bm{X}_1-\bm{X}_0\bm{W}\parallel$，再加上一个惩罚项}，该惩罚项应该是$ \bm{X}_1 $与通过权重$ \bm{W} $组合的相应控制组的值之间距离的增函数。或者就像前面说的，缩小待选池。

\textbf{\color{RoyalBlue}如何选择$ \bm{V} $}？任何$ \bm{V} $都不会影响我们推断程序有效，但是一个好的$ \bm{V} $还是可以减少合成控制估计量的均方误。有几个办法：
\begin{itemize}
	\item 利用$ \bm{X}_0,\bm{X}_1 $的预测能力的主观评价；
	\item 在干预前时期，合成控制地区（即未干预地区）可以逼近干预地区结果变量。如根据Abadie and Gardeazabal(2003)，在干预前时期，选一个正定对角矩阵使得因变量的均方误最小；
	\begin{equation}\label{6}
	\arg\min _{V\in\mathcal{V}}(\bm{Y}_1-\bm{Y}_0\bm{W}^*(\bm{V}))'(\bm{Y}_1-\bm{Y}_0\bm{W}^*(\bm{V}))
	\end{equation}
	其中$ \bm{Y}_0,\bm{Y}_1 $是干预前时期某个子集的控制组和处理组结果变量。这样，通过\eqref{5}式，任选一个$ \bm{V} $，就有一个最优的$ \bm{W} $，也即形成函数关系$ \bm{W}^*(\bm{V}) $，再通过\eqref{6}式就可以考察这个$ \bm{W} $是不是最优化了\eqref{6}式，如果不是，继续选择$ \bm{V} $，直到找到。
	\item 若样本中干预前时期足够长，则可将样本划分为训练样本与预测样本，对于任意的$ \bm{V} $，可以从训练样本中得到$ \bm{W}^*(\bm{V}) $，则可以通过在预测样本中计算，选一个由于$ \bm{W}^*(\bm{V}) $所导致的最小的预测均方误的$ \bm{V} $。
\end{itemize}

\subsubsection{推断}
基于回归的比较案例研究通常报告报告的标准误仅仅反映加总数据的不确定性（如Card,1990; Card and Krueger,1994）,另外，大样本推断技术也不适合比较组中的个体比较小的情况。本文提出的推断技术类似于\textbf{\color{RoyalBlue}置换检验(permutation test)}，无论数据是加总还是个体，也不要求比较组的数目很多。

检验统计量的计算是将干预组和非干预组进行了随机安排，这就可以让我们评价一个真实干预组的效应和随机选择一个组作为干预组的效应是否相对更大。一般来说，真实干预的效应应该相对没有干预的地区估计的效应的分布要相对大。若干预组较少，则可以利用长时间序列进行安慰剂检验（Mullainathan,2004），可在干预前选任意个时间作为伪干预时间。

\subsection{实例：估计California99项目}
\subsubsection{前期说明}
\begin{table}[H]
	\begin{tabular}{lp{30em}}
		背景：& 99项目是反吸烟，1988年在加利福利亚州开始立法，1989年1月生效；\\
		数据：&1970-2000年州水平的面板数据，干预前年数是19年。现在，我们来重现加利福利亚的反事实构造，有四个州因为滞后也采取了香烟控制，所以这四个州从待选池中删掉。同时，那些把香烟税提高了50美分的州也删掉。因此，待选池有38个美国州;\\	
		结果变量：&州水平的人均香烟消费；\\
		$ \bm{X}_0,\bm{X}_1 $中的解释变量：&平均香烟零售价、人均州个人收入对数、15-24岁人口比例、人均啤酒消费、三个滞后因变量（1975、1980和1988）；（当加上其他解释变量，如失业、收入不平等、贫穷、福利转移、犯罪率、香烟税、人口密度等也不影响结果）
	\end{tabular}
\end{table}
\subsubsection{基本结果}
\begin{figure}[H]
	\includegraphics[height=8cm]{f1.png}\hspace{2em}
	\includegraphics[height=8cm]{f3.png}
\end{figure}

左图实线是加利福利亚州人均香烟消费的时间序列，虚线是按照人口加权平均的其他控制组的的人均香烟消费。右图的虚线是合成控制后的对应时间序列。
\begin{figure}[H]
	\includegraphics[scale=0.7]{f2.png}
\end{figure}
上表第一列是加利福利亚州的解释变量均值，第二列是合成控制的解释变量均值，最后一列是人口加权的解释变量均值。
\subsubsection{推断}
迭代构造一系列安慰剂检验：即每次迭代就是从38个控制州中拿一个出来当做干预组，把加利福利亚州放入控制组。这样构造出来的效应的分布如下图，
\begin{figure}[H]
	\includegraphics[height=8cm]{f4.png}
\end{figure}


考虑到若干预前的拟合较差，那么干预后出现的巨大差距也不应看作是干预效果，或许本身就是拟合差的结果。因此，\textbf{\color{RoyalBlue}要把拟合比较差的州删掉，下面图中的左上、右上、左下三个子图分别是删掉20倍、5倍和2倍于加利福亚州均方误的图}。最后一个图计算的是干预后均方误除以干预前均方误的比例，可以看到加利福利亚州的最大。
\begin{figure}[H]
	\includegraphics[height=8cm]{f5.png}\hspace{2em}
	\includegraphics[height=8cm]{f6.png}\hspace{2em}
	\includegraphics[height=8cm]{f7.png}\hspace{2em}
	\includegraphics[height=8cm]{f8.png}	
\end{figure}

\subsection{多个干预组：Acemoglu et. al.（2016）}
前面看到处理组只有一个，控制组可以有多个。如果处理组也有多个，怎么估计呢？Acemoglu et. al.(2016)有一个简单的办法，令合成控制构造的结果变量为，
\[ \widehat{Y}_{it}=\sum_{j\in\text{Control group}}w_j^iY_{jt}  \]

其中$ w_j^i $是针对处理组个体$ i $的$ j $个控制组个体对应的权重。此时，干预效果为，
\[ \widehat \phi(\tau,k)=\frac{\sum_{i\in\text{处理组}}\frac{\sum_{t\in\text{处理后时期}}Y_{it}-\widehat Y_{it}}{\hat \sigma_i}}{\sum_{i\in\text{处理组}}\frac{1}{\hat \sigma_i}} \]

其中，
\[ \widehat{\sigma}_i=\sqrt{\frac{\sum_{t\in \text{处理前时期}}(Y_{it}-\widehat Y_{it})^2}{T}}\]


这种方法的本质就是一种加权，更好的拟合有更大的权重。

\subsection{回归合成法}
Hsiao et al. (2012) 提出了一种新的合成方法，在理论证明上，其另起灶炉，但与合成控制法最大的区别在于，基于事前时期用控制组个体合成干预组个体时，其权重可以为负，可以有截距项。

在事前时期远远大于控制组个体数时，可以选所有控制组个体作为潜在合成控制。否则，估计精度会降低，Hsiao et al. (2012)提出一种两步法解决该问题（穷举法）：
\begin{itemize}
	\item 第一步，在只有一个控制组个体时，可以估计$ C_N^1=N $个模型，根据$ R^2 $或似然值挑一个最好的，记为$ M(1)^* $。只有两个控制组个体可以估计$ C_N^2=2!/[N!(N-2)!] $个模型，也可以挑一个最好的，记为$ M(2)^* $，依次类推。
	\item 第二步，对于第一步得到的$ M(1)^*,M(2)^*,\cdots,M(N)^* $，根据AIC等准则再挑一个最好的即可。
\end{itemize}

\subsection{R包：合成控制Synth包}
文章作者也写了一个R包Synth，可以用来进行合成控制。一般其命令过程为：
\begin{itemize}
	\item 调用\verb|dataprep()|形成$ \bm{X}_0,\bm{X}_1,\bm{Z}_0,\bm{Z}_1 $（\textbf{\color{RoyalBlue}注意包中的$ \bm{X}_0,\bm{X}_1,\bm{Z}_0,\bm{Z}_1 $与论文中相应字母的含义不同}，包中分别指的是控制组解释变量，处理组解释变量，控制组因变量和处理组因变量）；
	\item 调用\verb|synth()|进行合成估计；
	\item 调用\verb|synth.tab(),synth.plot(),path.plot()|查看估计结果；
\end{itemize}


\subsubsection{命令的相关参数}
\verb|dataprep()|包含如下参数：
\begin{table}[H]
	\begin{tabular}{lp{30em}}
		\verb|foo|:& 原始的数据框；\\
		\verb|predictors|:& 自变量名，字符型；也就是下图中的\verb|X1,X2,X3|。\\
		\verb|predictors.op|:& 如何将不同时期的$ X_{it} $整合成一个$ \bar X_i $，类似于将$ Y $换成$ Z $时的\eqref{4}式，$ X $也需要一个这样的转换。譬如均值\verb|mean|。\\
		\verb|special.predictors|:& 类似文献中的滞后因变量设置或者其他特殊自变量的设置（比如处理方法\verb|predictors.op|不是均值啊之类的），这是一个列表，分别包含变量名、区间段和方法（类似\verb|predictors.op|）。\\
		\verb|dependent|:& 因变量名，字符型；即下图中的\verb|Y|。\\
		\verb|unit.variable|:& 面板中标识个体的id，即下图中的\verb|unit.num|，注意该变量应该是数值型。\\
		\verb|time.variable|:&面板中标识时间的id，即下图中的\verb|year|，注意该变量应该是数值型。\\
		\verb|treatment.identifier|:&即下图中变量\verb|unit.num|中的\textbf{\color{RoyalBlue}7}。\\
		\verb|controls.identifier|:& 即下图中变量\verb|unit.num|中除\textbf{\color{RoyalBlue}7}以外的如2等。\\
		\verb|time.predictors.prior|:& 数值型向量，标识处理前时期。\\
		\verb|time.optimize.ssr|:& 标识时间的数值型向量。这段时间的处理组因变量和估计的合成控制组因变量其均方误要最小。\\
		\verb|time.plot|:& 数值向量，标识要画出哪些时间的图\\
		\verb|unit.names.variable|:&  面板中标识个体的id，与\verb|unit.variable|不同的是它的值是字符而非数值型，即下图中的\verb|name|。
	\end{tabular}
\end{table}

数据框例子：
\begin{figure}[H]
	\centering
	\includegraphics[height=10cm]{f9.png}
\end{figure}

该函数输出$ \bm{X}_0,\bm{X}_1,\bm{Z}_0,\bm{Z}_1 $等。

\vspace{5em}
\verb|synth()|包含如下参数：
\begin{table}[H]
	\begin{tabular}{lp{35em}}
		\verb|data.prep.obj|:&	就是\verb|dataprep()|命令的输出；\\
		\verb|X1,X0,Z1,Z0|:& 当有了\verb|data.prep.obj|参数时，这些参数就无需输入；\\
		\verb|custom.v|:&	该矩阵可以由你指定，也可以不输入，由命令自行优化运算，注意，自行指定了\textbf{V}，实际上就已指定了\textbf{W}；\\
		\verb|optimxmethod|:&	指定使用何种优化算法，是拟牛顿啊还是其他之类的；\\
		\verb|genoud|:& 逻辑标识。为真，则意味着进行两段优化，第一阶段通过进化算法找一个大约的值，第二阶段在这个值附近进行优化算法找到最优值。一般会耗费更多的时间。	\\
		\verb|quadopt|:& 设定二次规划规则，缺省的是\verb|ipop|；\\
		\verb|Margin.ipop|:& 二次规划时离约束有多近；\\
		\verb|Sigf.ipop|:& 二次规划的精度，几位有效数字；\\	
		\verb|Bound.ipop|:&二次	规划的修剪边界；\\
		\verb|verbose|:&是否展现中间过程；
	\end{tabular}
\end{table}

\subsection{R包：回归合成pampe包}
譬如Hsiao et al. (2012)估计香港回归后的经济增长效应，那么数据如下排列：
\begin{figure}[H]
\includegraphics[scale=0.6]{pamp_1.png}
\end{figure}

pampe函数的参数：
\begin{table}[H]
\begin{tabular}{lp{40em}}
time.pretr	& 干预前时期（不包括干预当期），如1：10\\
time.tr	& 干预后时期，如11：20\\
treated & 处理个体，即列名或数值型第几列\\
controls & 同上，但是控制个体.缺省的是设置了处理个体，剩下的都是控制个体\\
data& 可以是数据框，格式见上图，行是时间，列是待选控制组个体\\
nbest & 选模型时根据$ R^2 $每一步选几个最好的模型出来。可以选2个，3个等。\\
nvmax & 选多少个控制组进来检查模型拟合优度。缺省是全部	，如果数目太多会抛错。\\
select & 第二步用什么准则选模型，如\verb|AIC|等\\
placebos & 安慰剂检验，三个可选："controls", "time", or "Both"。"controls" 重新安排干预给控制组，计算因果效应。"time"在真实干预之前，重新安排干预时间。
\end{tabular}
\end{table}


\section{RDD：断点回归设计}
有一个结果变量$ Y $，还有一个参考变量$ X $，若$ X\ge x_0 $，则一定被干预，即$ D_i=1(X\ge x_0) $，这被称为\textbf{\color{RoyalBlue}精确断点}。如果若$ X\ge x_0 $，只是被干预的可能性增大，不一定被干预，这称为\textbf{\color{RoyalBlue}模糊断点}。如凭学生成绩拿奖学金，若上了线一定能拿奖学金，就是精确断点。上了线，还有其他考量，只是上线以后拿奖学金概率增大，那么就是模糊断点。

\subsection{假设条件}
\begin{itemize}
	\item 断点假设：若$ p^+=\lim\limits_{x\rightarrow x_0^+}E[D_i|X_i=x],p^-=\lim\limits_{x\rightarrow x_0^-}E[D_i|X_i=x] $存在，且
	\[ p^+\ne p^- \]
	通俗地说，个体分配概率在临界值左右有跳跃。
	\item 连续性假设：令$ E[Y_{0i}|X_i=x],E[Y_{1i}|X_i=x] $是$ x $的函数，并且在$ x_0 $处是连续的，即，
	\[ \lim\limits_{\varepsilon\rightarrow 0}E[Y_{ji}|X_i=x_0+\varepsilon]=E[Y_{ji}|X_i=x_0-\varepsilon],\qquad j=0,1 \]
	\item 局部随机化假设：在断点附近近似于完全随机化实验，即
	\[ (Y_{0i},Y_{1i})\perp D_i|X_i\in \delta(x_0) \]
	其中$ \delta(x_0)=(x_0-\delta,x_0+\delta)$为$ x_0 $的邻域，$ \delta>0 $为任意小的正数。即不能精确控制或操纵参考变量使之超过临界值。比如高考上一本线，对于在一本线附近的同学，其上线存在着很大运气成分。	
\end{itemize}

根据上述三个假设，\textbf{\color{RoyalBlue}因果效应$ \tau_i=Y_{1i}-Y_{0i} $的估计为(Hahn et., 2001)}，
\[ E[\tau_i|X_i=x_0]=\frac{\mu^+-\mu^-}{p^+-p^-}=\mu^+-\mu^- \]

其中，$ \mu(x)=E[Y_i|X_i=x],Y_i=Y_{0i}+\tau_iD_i,\mu^+=\lim\limits_{x\rightarrow x_0^+}\mu(x), \mu^-=\lim\limits_{x\rightarrow x_0^-}\mu(x)$

但在模糊断点下，局部随机化假设往往不能满足，需要其他的替代假设：定义$ D_i(x) $为参考变量为$ x $时的个体干预状态，$ D_{1i}(x) $为参考变量在断点右侧的个体参与状态，$ D_{0i}(x) $为参考变量在断点左侧时的个体参与状态，即，
\[ D_i=\begin{cases}
D_{1i}\qquad x\ge x_0\\
D_{0i}\qquad x< x_0
\end{cases} \]

\begin{itemize}
	\item 独立性假设：$ Y_{0i},Y_{1i},D_{1i},D_{0i} $在断点附近独立于参考变量$ X_i $，即
	\[ (Y_{1i},Y_{0i}),D_{1i},D_{0i}\perp X_i,X_i\in\delta(x_0) \]
	\item 单调性假设：断点对个体的影响方向是一致的，对于任意的$ x\in \delta(x) $，若假设正向单调性成立，有，
	\[ D_{1i}(x)\ge D_{0i}(x) \]
	即个体处在右侧时总是比处在左侧时更要被干预。
\end{itemize}

用这两个假设替代局部随机化假设，有（Hahn et., 2001），
\begin{equation}\label{RDD_1}
 E[\tau_i|D_{1i}(x)>D_{0i}(x)]=\frac{\mu^+-\mu^-}{p^+-p^-} 
\end{equation}


\subsection{RDD的实施}
\begin{itemize}
	\item 第一步，画图，检查变量的连续性。
	\item 第二步，估计。
	\item 第三步，稳健性检验。
\end{itemize}
\subsubsection{第一步：画图}
为何要画图？干预分配概率与结果变量在临界点会有跳跃，而其他协变量没有跳跃。因此，
\begin{itemize}
	\item 可以绘制参考变量与结果变量的关系图，观测结果变量在断点处是否有跳跃。也可以看看它在其他位置是否也是这样。直接绘图，可能会包含很多噪音，可类似于核密度估计平滑后，绘制。
	\item 可以绘制其他协变量与参考变量的关系图，观其是否连续。
	\item 绘制参考变量的分布图，观其分布是否在断点处连续，若不连续，则说明个体能精确控制参考变量。光看图不精确，McCrary (2008) 还提出了一个检验方法。
\end{itemize}

\subsubsection{第二步：估计}
多种估计方法：\textbf{\color{RoyalBlue}边界非参数回归、局部线性回归、局部多项式回归}。因非参回归收敛速度慢，Imbens and Lemieux (2008)建议采用局部线性回归。
\paragraph{边界非参数估计}
譬如现有一组数据$ (X_i,Y_i) $，此时给你一个值$ x $，对应它的$ y $是多少呢？最简单的思路就是，在$ x $左右某个区间$(x-h,x+h)$内观测值的$ Y $点取平均即可，即，
\[ \bar y=\frac{1}{N}\sum_{i=1}^{N}Y_i\cdot \bm{1}(x\in [x-h,x+h]) \]

进一步地，如果我们想给予接近$ x_k $更大的权重，那么就不能像上式一样等权重求和，而是加权求和，选择什么样的权重就是选择什么样的核，
\[ \bar y=\frac{Y_1\cdot K(\frac{X_1-x}{h})}{\sum_{i=1}^{N}K(\frac{X_1-x}{h})}+\cdots+ \frac{Y_N\cdot K(\frac{X_N-x}{h})}{\sum_{i=1}^{N}K(\frac{X_N-x}{h})}=\sum_{i=1}^{N}\frac{Y_i\cdot K(\frac{X_i-x}{h})}{\sum_{i=1}^{N}K(\frac{X_i-x}{h})}\]

这样，通过在断点左右两边选择非参数估计，然后求差来得到因果估计。

\paragraph{局部线性回归}边界非参数估计至少有两个缺点：
\begin{itemize}
	\item 局部随机化假设要求在断点附近很小的范围，而在此范围，边界非参数估计可能没有足够的样本点，这就会带来估计的偏误.
	\item 若参考变量$ X $对潜在结果也有影响，那么考虑$ X $后再估计会更好。
\end{itemize}。
考虑局部线性回归进行调整，
\begin{align*}
\min_{a_l,b_l} \sum_{i=1}^N[Y_i-a_l-b_l(X_i-x_0)]^2K\left(\frac{X_i-x_0}{h}\right)\cdot \bm{1}(X_i<x_0)\\
\min_{a_h,b_h} \sum_{i=1}^N[Y_i-a_h-b_h(X_i-x_0)]^2K\left(\frac{X_i-x_0}{h}\right)\cdot \bm{1}(X_i\ge x_0)
\end{align*}

上述两个方程对左右两边均值的估计为，
\begin{align*}
\hat \mu^-(x_0)=\hat a_l+\hat b_l(x_0-x_0)=\hat a_l\\
\hat \mu^+(x_0)=\hat a_h+\hat b_h(x_0-x_0)=\hat a_h
\end{align*}

因此，两个回归的截距之差就是因果效应。可以看到局部线性调整和边界非参数估计的差别在于对谁加权，而且局部线性还有一个最小化选择系数的过程。\textbf{\color{RoyalBlue}如果范围较窄，线性回归可以对潜在结果与参考变量间的非线性形式进行较好的近似。}

实际上，这个因果效应也可以通过如下方式直接得到（Imbens and Lemieux, 2008），
\[\min_{a,b,\tau,\gamma}\sum_{i=1}^N[Y_i-a-b(X_i-x_0)-\tau D_i-\gamma D_i(X_i-x_0)]^2K\left(\frac{X_i-x_0}{h}\right)\cdot \bm{1}(x_0-h\le X_i\le x_0+h)\]

\paragraph{局部多项式回归}当样本量实在太少，可能需要加宽$ h $，此时可以考虑局部多项式估计以捕捉可能的非线性关系，此时就在回归中添加高次项即可。

\paragraph{模糊断点回归的估计}
注意到\eqref{RDD_1}式是一个比值，清晰断点可以求得分子，实际上分母可以通过令因变量为$ D_i $(表示处理状态，此时不是临界值的确定性函数，超越了断点可能处理也可能不处理)作同样的上述三种算法来得到，譬如，若使用局部线性算法，

那么可以使用，
\[ \min_{a_Y,b_Y,\tau_Y,\gamma_Y}\sum_{i=1}^NK\left(\frac{X_i-x_0}{h}\right)\cdot [Y_i-a_Y-b_Y(X_i-x_0)-\tau_YT_i-\gamma_YT_i(X_i-x_0)]^2 \]
以得到$ \tau_Y $，其中$ T_i=\bm{1}(X_i\ge x_0) $。使用，
\[ \min_{a_Y,b_Y,\tau_Y,\gamma_Y}\sum_{i=1}^NK\left(\frac{X_i-x_0}{h}\right)\cdot [D_i-a_Y-b_Y(X_i-x_0)-\tau_YT_i-\gamma_YT_i(X_i-x_0)]^2 \]
以得到$ \tau_D $，那么模糊断点估计量为，
\[ \color{RoyalBlue}\tau^{FGD}=\frac{\hat \tau_Y}{\hat \tau_D} \]
\paragraph{如何选择带宽？}
涉及到非参数估计或者局部线性（非线性）回归（\verb|DifferentBinwith.R|），就一定会有最优带宽的选择。如何选呢？
\begin{enumerate}
	\item 只针对临界点左右两边一定范围之内的样本进行运算、搜索；
	\item 对某个带宽$ h $，总可以得到左侧（右侧）回归参数$ a_l,b_l(a_h,b_h) $的估计，利用该估计，针对左侧（右侧）的任意样本点，总有一个预测值$ \hat\mu(X_i) $，使得该预测值与原始值$ Y_i $之间的距离最小的$ h $就是最好的带宽。即，
	\[ h^{opt}=\arg\min_h \frac{1}{N}\sum_{i=1}^N[Y_i-\hat\mu(X_i)]^2 \]
	\item 对于模糊断点，需要对分母也进行一次带宽选择，$ h^{opt}_D=\arg\min_h \frac{1}{N}\sum_{i=1}^N[D_i-\hat\mu(X_i)]^2 $。Imbens and Lemieus (2008) 建议分子分母使用相同带宽，用最小的带宽即可。
	\item 带宽的选择，Imbens and Kalyanaraman (2012)（IK带宽，代码中有）， Calonico et al. (2014)（CCT带宽）, Calonico et al. (2017)又进行了改进。
\end{enumerate}

\subsubsection{第三步，稳健性检验}
\begin{itemize}
	\item 协变量连续性检验，其实就是绘制协变量与参考变量的关系图；
	\item 伪断点检验：另外选个附件的点，看是不是有因果效应。譬如断点左右两侧的中点位置。
	\item 带宽的敏感性检验。不同带宽对RDD估计量是否产生重大影响。
\end{itemize}

\verb|rdrobust(y = vote, x = margin)|
\subsection{实践}
主要调用\verb|rddapp|包，估计函数为\verb|rd_est|，其相关参数如下：
\begin{table}[H]
\begin{tabular}{lp{40em}}
formula	& 简单的精确断点，直接写为\verb|y ~ x|，若还有两个协变量就写为\verb"y ~ x | c1 + c2"。模糊断点写为，\verb|y ~ x + z|，\verb|x|是参考变量，\verb|z|是内生干预变量（即$ D_i $）。\\
data & 	数据框\\
subset	& 一个向量以选定子集\\
cutpoint & 断点是多少，忽略的话，断点是0\\
bw & 带宽，数值型。给定一个带宽，它会按它的一半，它的两倍，它本身三种带宽输入估计结果。缺省的是“IK12”方法选择的带宽\\
kernel	& 局部线性回归中选择的核。包括缺省的三角核"triangular"， 矩形核"rectangular",以及其他核"epanechnikov", "quartic", "triweight", "tricube", "gaussian" and "cosine"\\
se.type & 标准误的计算方法，包括"HC3", "const", "HC", "HC0", "HC1", "HC2", "HC4", "HC4m", "HC5"\\
cluster & 误差的集聚方式，影响标准误计算\\
t.design & 干预的方式：参考变量大于断点则被干预，选"g"； 参考变量大于或等于断点则被干预，选"geq"；参考变量小于断点则被干预，选 "l" ；参考变量小于或等于断点则被干预，选 "leq" 
\end{tabular}
\end{table}

\verb|rdrobust|包绘断点图很不错。
%\subsection{Malmquist生产指数: \texttt{productivity}}
%\section{IV}

\include{TS}


\end{document}